/**********************************************************
 *
 *  Copyright (c) 2003  SeikoEpson Inc
 *  All Rights Reserved
 *
 *  File name : f_fmod.c
 *  Function  :
 *        f_fmod function returns a remainder when deviding 
 *        the value.
 *        This file is copied from math.lib of 
 *                  CC33 tool(CC33v40.exe + CCPA45.exe, 
 *                  math of ansilib33v25 )
 *  original file's Revision  :
 *      2000/02/14    first release                 M.Igarashi
 *  Revision  :
 *      2003/04/09    IrumaSoft M.Takeishi   1.st design
 *
 **********************************************************/
#include <f_math.h>
#include <f_smcvals.h>
 
// ALGORITHM
// 1.|sfX|>|sfX2| fmod (sfX,sfX2)= sfX2 * modf( sfX / sfX2,&d)
//   |sfX2|>|sfX| fmod (sfX,sfX2)= sfX
//   |sfX|=|sfX2| fmod (sfX,sfX2)= 0
//

//
//  the architecture of float floating point
//
//   0 30         23 22                0  32 bit
//   -----------------------------------
//  | |   exponent  |  fraction         |
//   -----------------------------------
//
//  |               |                   |
//  |    8 bits     |    23 bits        |
//
//         bit    31         sign bit        (  1 bit  )
//             30 - 23      exponent part    (  8 bits )
//             22 -  0      fraction part    ( 23 bits )
//
//
//

//#define FOR_MORE_PRECISION
//   If you define FOR_MORE_PRECISION, you'll get more precision, but
//   the run time increases. But sometimes f_fmod returns a very bad error 
//   value in defining FOR_MORE_PRECISION.
//

float f_fmod (float sfX, float sfX2){

	long lX,lX2;
	unsigned long ulx2;
	float sfTemp,sfRet,sfMod;

#ifdef FOR_MORE_PRECISION
	FLT_LNG sfTemp2;
#endif

	F_GETW(lX,sfX);			// get low
	F_GETW(lX2,sfX2);		// get low


	ulx2=lX2&0x7fffffff;	// mask sign
	if(ulx2==0x00000000 ){		// sfX2 ==0
	    errno=EDOM;                         // 33: domain error
	    // performance tunning
	    //sfRet= f_NAN._F; 
		//return sfRet;
		return f_NAN._F;
	}

	/* first arg , NaN, Inf check */
	ulx2= 0x7f800000;	// mask sign
	if ( ( lX & ulx2 ) == ulx2 ) {

		ulx2 = 0x7fc00000;
		if ( ( lX | 0x80000000 ) == 0xff800000 || 
		    ( lX & ulx2 ) == ulx2 ) {
			/* +-Inf or NAN */
			errno = EDOM;
			sfRet = f_NAN._F;
			return sfRet;
		}
	}

	/* 2nd arg, NaN check */
	ulx2 = 0x7fc00000;
	if ( ( lX2 & ulx2 ) == ulx2 ) {
		errno = EDOM;
		sfRet = f_NAN._F;
		return sfRet;
	}


	if(((lX2&0x80000000) !=0 ? -sfX2:sfX2) < ((lX&0x80000000) !=0 ? -sfX:sfX)){		// |sfX|>|sfX2|


		sfTemp=sfX/sfX2;

		sfTemp=f_modf(sfTemp,&sfMod);

#ifdef FOR_MORE_PRECISION
		// get ignored value, and get its fraction part
		sfTemp2._F = sfX - ( ((float)sfMod) * sfX2 );
		sfTemp2._F -= sfTemp * sfX2;
		if ( sfTemp2.st._LL != 0x00000000 ) {
			sfTemp2._F /= sfX2;
			sfTemp2._F =f_modf(sfTemp2._F ,&sfMod);
		} else {
			sfTemp2.st._LL = 0x00000000;
		}
#endif

		sfRet=sfX2*sfTemp;
#ifdef FOR_MORE_PRECISION
		if ( sfTemp2.st._LL != 0x0 ) {
			sfTemp2._F *= sfX2;
			sfRet += sfTemp2._F;
		}
#endif


		return sfRet;
	}else{
		if(((lX2&0x80000000) !=0 ? -sfX2:sfX2) > ((lX&0x80000000) !=0 ? -sfX:sfX)){	// |dfX2|>|dfX|
			//sfRet=sfX;
			//return sfRet;
			return sfX;
		}else{				// dfX==dfX2
			//sfRet=0.0f;
			//return dfRet;
			return 0.0f;
		}
	}
}
